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SUMMARY 


Seismograms  were  synthesized  by  the  Gaussian  beam  method  and  ray  trajectories 
were  plotted  for  a  plane  wave  incident  on  a  3-D  structure  obtained  from  block  inversion 
of  teleseismic  travel  times  observed  at  a  local  array  in  central  California.  The 
amplitude  variations  for  the  given  profile  are  smaller  than  the  order  of  magnitude 
variations  typically  observed.  Since  amplitudes  are  sensitive  to  smaller  scale  velocity 
features  than  travel  times,  this  may  indicate  that  a  single  scale  block  model  may  have 
difficulty  in  satisfying  both  travel  time  and  amplitude. 

The  synthesis  of  teleseismic  body  waves  by  either  the  Gaussian  beam.  .Maslov,  or 
Kirchhoff  techniques  entails  the  integration  of  dynamic  and  kinematic  ray  tracing 
equations  over  parts  of  the  model  having  3-D  as  well  as  1-D  velocity  variations.  The 
longest  integration  intervals  are  over  the  1-D  or  radially  symmetric  part  of  the  model. 

A  procedure  consisting  of  the  multiplication  of  propagator  matrices  can  be  used  to 
smoothly  patch  analytic  solutions  of  the  dynamic  equations  obtained  over  the  1-D  or 
radially  symmetric  part  of  the  model  into  numerical  solutions  obtained  over  the  parts 
of  the  model  having  3-D  velocity  variations. 
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1.  INTRODUCTION 

Observations  of  teleseismic  P  waves  at  large  aperture  arrays  have  found  amplitude 
fluctuations  across  the  array  as  large  as  those  observed  over  the  world  wide  networks. 
The  amplitude  fluctuations  and  their  correlation  with  travel  time  variations  are 
consistent  with  the  focussing/defocussing  effects  of  three-dimensional  velocity 
structure  beneath  the  array  rather  than  the  effects  of  intrinsic  attenuation.  In  this 
semi-annual  report,  the  Gaussian  beam  method  is  used  to  predict  the  possible 
amplitude  variations  at  a  large  aperture  array  due  to  a  known  3-D  structure  The 
known  3-D  structure  is  determined  from  a  block-inversion  of  travel  times  of  the  type 
described  by  Aki  et  al.  (1976).  For  body  waves  on  the  order  of  1  Hz.  dominant 
frequency,  the  scale  length  and  intensity  variations  of  the  models  commonly  derived 
from  such  inversions  are  well  within  the  bounds  in  which  asymptotically  approximate 
methods,  such  as  the  Gaussian  beam  method,  are  valid. 

A  goal  of  these  experiments  is  to  determine  whether  knowledge  of  the  3-D 
structure  surrounding  a  given  nuclear  test  site  can  be  used  to  reduce  the  uncertainty 
in  a  yield  estimate  by  providing  an  amplitude  correction  that  is  a  function  of  source 
location  within  the  test  site.  Although  the  spatial  resolution  of  3-D  block  inversions  is 
much  less  than  that  obtainable  from  on-site  reflection  experiments  and  well-logs,  it  is  a 
less  intrusive  experiment.  Such  inversions  are  always  possible  if  the  events  within  a 
test  site  are  well  timed  and  located.  Thus  it  is  important  (l)  to  determine  whether  the 
spatial  resolution  of  3-D  block  inversions  is  sufficient  to  explain  any  significant 
amplitude  variations  in  teleseismic  P  waves  due  to  variations  of  event  locations  within  a 
lest  site  and  (2)  to  determine  the  relation  between  characteristic  scale  lengths  and 
intensities  of  velocity  fluctuations  of  these  models  with  predicted  amplitude 
fluctuations. 

The  preliminary  calculations  shown  in  this  report  have  approximated  the  effects  of 
velocity  variations  beneath  a  test  site  by  the  reciprocal  problem  of  a  plane  wave  front 
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incident  on  a^hree-dimensionally  varying  structure.  In  order  to  more  precisely 
calculate  the  teleseismic  problem,  it  is  necessary  to  determine  the  geometric  spreading 
and  Gaussian  beam  amplitudes  over  long  distances  through  parts  of  the  model  that  are 
likely  to  have  much  weaker  velocity  fluctuations.  These  parts  of  the  model  can  be 
represented  by  a  1-D  or  radially  symmetric  earth.  In  these  parts  of  the  model,  dynamic 
ray  tracing  and  quantities  needed  for  Gaussian  beam  or  Maslov  synthesis  can  be  rapidly 
calculated  using  a  propagator  matrix  technique.  Section  4  of  this  report  gives  analytic 
forms  for  the  elements  of  this  propagator  in  a  1-D  earth  model. 

3  THE  GA I  'SSI AN  BEAM  METHOD  OF  SEISMOGRAM  SYNTHESIS 

The  Gaussian  beam  method  is  an  extension  of  the  ray  method  for  the  computation 
of  seismograms  in  a  smoothly  varying  heterogeneous  medium  {see  f^erven^  et  al.,  19B2). 
In  comparison  to  other  methods,  the  Gaussian  beam  method  is  fast,  gives  finite 
amplitudes  at  caustics,  and  requires  no  explicit  two-point  ray  tracing.  The  method  uses 
rays  as  a  framework  upon  which  the  wavefield  is  built.  The  Gaussian  beams  are 
propagated  along  each  ray  and  then  superposed  at  the  receiver.  The  beams  are 
weighted  in  the  superposition  in  order  to  satisfy  a  given  source  condition. 

The  individual  beams  in  the  superposition  can  be  manipulated  by  a  complex  beam 
parameter  which  changes  the  beam  curvature  and  beam  width  at  the  source.  This  is 
done  in  order  to  satisfy  different  objectives,  such  as  small  beam  width  at  the  source, 
receiver,  or  a’  a  localized  heterogeneity,  or  to  satisfy  a  given  curvature  condition.  It 
has  been  found  that  large  initial  beam  widths  are  more  suitable  for  the  decomposition 
of  a  point  source  into  Gaussian  beams  and  small  initial  beam  widths  are  more  suitable 
for  the  decomposition  of  a  plane  wave  into  Gaussian  beams  (see  Nowack  and  Aki,  1984). 
The  Gaussian  beam  method  reduces  to  a  plane  wave  expansion  for  large  planar  beams 
at  a  point  source  in  a  homogeneous  source  region.  Finite  beam  widths  have  the  effect 
of  localizing  the  beams  to  a  vicinity  of  each  central  ray,  and  reducing  numerical  end 
effects  in  the  superposition  (see  Madariaga  and  Fapadimitriou,  1905).  Finite  beam 
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widths  also  reduce  diffraction  effects  in  certain  cases,  such  as  the  generation  of  a  head 
wave  at  a  plane  interface.  The  critical  initial  beam  width  gives  the  beam  solution  most 
concentrated  about  each  central  ray  and  is  a  compromise  that  has  been  found  to  give 
stable  results.  Tests  at  a  simple  caustic  have  nonetheless  given  the  correct  behavior 
over  a  range  of  initial  beam  widths.  More  work  must  be  done,  particularly  in  1-D  and  2- 
D,  in  choosing  the  appropriate  beam  parameter  that  gives  the  best  results  for  a  given 
heterogeneous  medium. 

The  numerical  advantages  of  using  high  frequency  asymptotic  methods  are  most 
pronouru.’ed  in  three-dimensional  cases.  The  finite  difference  method  in  contrast  is 
slow  for  2-1)  wave  propagation  problems,  and  impractical  for  most  3-D  seismic  wave 
problems.  Still,  many  problems  require  approximate  3-D  solutions.  It  is  thus  of  interest 
to  develop  asymptotic  methods  such  as  the  Gaussian  beam  method  for  t  he  3-D  case 

3  RES  I  il  TS  WITH  A  3-D  MODEL 

An  example  3-D  structure  can  be  obtained  from  a  model  determined  from  a  travel 
time  inversion  parameterized  by  blocks  or  by  any  other  method.  Here,  the  results  from 
a  teleseismic  travel  time  inversion  study  for  central  California  by  Zandt  (19B1)  are 
used.  The  model  has  four  layers  from  0.0  to  90.0  km  in  depth.  The  horizontal  block  size 
i.s  10.0  km  in  the  top  layer  and  20.0  to  25.0  km  in  the  lower  layers.  Average  velocity 
variations  in  the  derived  model  are  between  4.0  to  B.O  percent  in  the  top  layer  and  2.0 
to  4.0  percent  in  the  lower  layers.  Comer  and  Aki  (19B2)  investigated  this  structure 
using  3-D  ray  tracing  in  order  to  study  the  effect  of  ray  bending  on  estimated  take-off 
angles  for  shallow  events.  Figure  la  shows  an  example  of  down  going  rays  from  a 
shallow  event  in  the  central  California  structure.  The  3-D  ray  tracer  used  is  a 
modification  of  a  subroutine  written  by  R.  Comer.  The  velocity  model  is  interpolated 
using  a  3-D  splines  under  tension  package  by  A.  Cline.  Figure  lb  shows  a  top  view  of  the 
ray  paths  in  which  ray  bending  can  clearly  be  seen.  The  rays  associated  with  a 
teleseismic  plane  wave  incident  from  below  is  shown  in  Figure  2a.  A  more  dense  set  of 
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rays  is  shown  in  Figure  2b.  Again,  ray  bending  can  be  seen  in  the  ray  diagram.s. 

To  compute  synthetic  seismograms,  the  paraxial  ray  equations  must  be  solved 
along  each  ray.  These  equations  are  a  linearized  set  of  equations  which  can  be  used  to 
obtain  information  in  the  vicinity  of  a  ray  that  has  already  been  traced  using  the 
kinematic  ray  equations.  The  paraxial  (dynamic)  ray  equations  are  useful  in  a  variety 
of  ways  including  the  computation  of  geometric  spreading  and  amplitude  in  the 
standard  ray  method,  the  calculation  of  local  wave  front  curvature  for  phase 
interpolation,  and  the  computation  of  polarization  vectors,  etc.  (see  Cerven;^,  19B5a). 
Complex  solutions  of  these  equations  are  used  in  the  Gaussian  beam  method.  Thus, 
once  the  kinematic  and  paraxial  ray  equations  are  solved  in  the  vicinity  of  the  receiver, 
this  information  can  be  input  directly  into  a  Gaussian  beam  synthesis  program. 

As  an  example,  Gaussian  beam  synthetics  are  computed  for  the  teleseismic  plane 
wave  incident  from  below  in  the  central  California  model.  The  endpoints  of  the  rays  are 
shown  in  Figure  3a.  A  2.0  Hz  Gabor  wavelet  is  used  as  the  source  time  function,  and  a 
critical  initial  beam  width  is  used.  A  correlation  between  lower  amplitudes  and  earlier 
arrivals  can  be  seen  on  the  synthetics.  The  lower  amplitude  can  be  identified  on  Figure 
3a  near  the  point  A.  where  the  rays  are  being  pulled  apart.  In  Figure  3b,  the  amplitude 
varies  by  a  factor  of  about  2  along  the  profile.  Observations  of  amplitude  at  the 
Montana  LASA  by  Aki  (1973)  showed  variations  in  the  amplitude  by  a  factor  of  10.  This 
difference  in  the  amplitude  variation  may  be  related  to  the  block  sizes  used  to  satisfy 
the  travel  time  data  for  the  inverted  block  model.  Since  travel  time  is  related  to  the 
velocity  field  and"  amplitude  is  related  to  the  second  derivative  of  the  velocity  field, 
different  heterogeneity  scales  may  be  required  to  satisfy  both  travel  time  and 
amplitude  variations. 


4.  THE  PROPAGATOR  MATRIX  FOR  DYNAMIC  RAY  TRACING 


Since  the  dynamic  (paraxial)  ray  tracing  system  can  be  written  in  the  linear  form 

=  SW.  it  permits  the  definition  of  fundamental  matrix  and  propagator  matrix 
as 

solutions.  Specifically,  Cerven^f  (1985a)  has  defined  a  4x4  matrix  n(s,So)  that  solves 
this  linear  system  by  four  linearly  independent  solutions,  with  s  denoting  a  point  along 
the  ray  path  where  the  integration  terminates  and  s,  a  reference  point  (t:.g.,  source 
point,  receiver  point)  where  the  integration  starts  and  initial  conditions  are  specified, 
n  has  the  property  that  n(s,  ,5^)  =  I ,  where  I  is  the  identity  matrix.  Also,  given  a  point 
s  between  s^  and  s  ,  n(s,So)  =  n(s,s  )  n(s  ,s<,)  .  This  property  can  be  used  to  accelerate 
two  point  ray  tracing  and  the  calculation  of  spreading  functions  and  dynamic  quantities 
needed  for  synthesis  of  body  waves  in  three-dirnensionally  varying  media  by  the 
Gaussian  beam  and  KirchhofT  techniques.  For  example,  if  the  three-dimensional 
variations  of  a  model  are  confined  only  to  the  portion  of  the  ray  path  between  s  and 
numerical  integration  of  the  linear  system  need  only  be  performed  between  Sg  and  .s  to 
construct  n(s  ,So)  .  Along  the  segment  of  the  ray  path  between  s  and  s,  traversing  a  1- 
D  or  homogeneous  portion  of  the  model,  analytic  solutions  exist  for  ri(s,.s  )  .  Thus 
n(.s,Sp'  an  be  simply  obtained  by  multiplying  the  analytic  and  numerically  obtained  FI 
matrices. 


The  n  matrix  can  be  defined  in  terms  of  x  2  sub  matrices  P  and  Q  as  follows: 


n  = 


Q  Q 

.P  p| 


(i: 


where  the  bar  over  Q  and  P  distinguishes  plane  wave  solutions  from  point  source 
solutions  to  the  dynamic  ray  tracing  equations.  The  dynamic  ray  tracing  equations  are 
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The  subscripts  on  the  elements  of  V  denote  differentiation  of  velocity  v  with  respect  to 
the  ray  centered  co-ordinate  directions  qi  and  gg.  In  order  to  determine  the  elements 
of  the  II  matrix,  one  must  first  select  a  vector  basis  of  the  ray  centered  co-ordinate 
system  (  gi.g2,f  )  at  the  point  at  which  itial  conditions  are  specified,  liere,  given  a  fixed 
('artesian  co-ordinate  system  (x.  y,  z),  the  ray  will  be  assumed  to  be  confined  to  the 
(y,/)  plane  in  a  1-D  medium  and  the  e,  dir(H‘tion  will  be  chosen  to  point  in  the  y 
direction,  .\  vector  eg  is  chosen  such,  that  (  e,e2.(  )  forms  an  orthogonal,  rigtii-lianded 
vi.'clor  liasis,  v.here  t  is  the  ray  tangimt. 

With  this  definition  of  the  ray  centered  co-ordinate  system,  and  assuming  tnat 
velocity  varies  only  as  a  function  of  z  (1-D  medium),  it  is  possible  to  der  ive  analy  tic 
expressions  for  the  elements  of  0.  Tin;  solutions  for  the  spreading  of  a  line  source  m  a 
vertically  {/)  varying  medium  have  been  given  by  Madariaga  (1934).  These  arc 
appropriate  for  the  22  components  of  the  Q  and  P  matrices.  For  the  spreading  of  a 
point  source  in  three  dimensions,  solutions  for  the  1 1  component  s  must  be  added  tc; 
these.  The  1 1  components  are  quite  simple  in  form.  This  is  because  by  the  choice'  of 
the  ray  centered  co-ordinates  in  the  1-D  medium,  n,)  =  0.  Also  the  off  diagonal 
element  s  of  the  Qand  P  mat  rices  are  all  zero  because  the  velocity  v  is  assumed  not  to 
vary  in  the  x  and  y  directions.  In  summary,  the  16  components  of  the  fl  matrix  in  a 
vertically  varying  medium  are 

4'n  =  =  1 
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Using  this  definition  of  the  propagator  FI,  it  can  be  verified  that  for  any  point  s 
along  the  ray  path  at  which  the  velocity  and  its  first  spatial  derivatives  are  continuous. 


X(s)  =  n(s.s)-n(s>o)-X(So)  (i 

where  X{s„)  is  a  2  x  4  matrix  of  initial  conditions  on  P  and  Q  and  X(s)  is  a  2  x  4  matrix 

formed  from  a  linear  combination  of  plane  wave  and  point  source  solutions  for  P  and  Q 

0  0 
0  0 

For  example,  one  can  show  that  given  point  source  initial  conditions  at  s,,  X(so)  =  j  q 


Q 

that  X{s)  =  pj ,  where  Q,  P  are  point  source  solutions  at  s  ,  and  also  that 

n(.v..Vo)  =  n(s,s  )n(s,So)  .  AIso,  using  this  definition  of  Flit  can  be  verified  that  the 
determinant  of  FI  is  constant  along  the  ray  path  and  that  det  FI  =  1 . 

If  the  velocity  or  its  first  spatial  derivative  are  discontinous  at  s’  ,  then  the  jump 
conditions  defined  by  Cerveny'  (1985a)  must  first  be  applied  to  the  elements  of  n(s  .So) 
before  the  multiplication  of  eq.  (3)  is  performed.  For  the  co-ordinate  system  chosen 
and  for  velocity  depending  only  on  z,  these  jump  conditions  are 


cusu 


cosrf  cosd*^ 


®(cos(5  V*  -  cogiS  v~)  Q^2 


(4b) 


P2Z 


cos6 


-^22 


+ 


where  the  +  sign  refers  to  quantities  on  the  incident  side  of  the  boundary  and  the  -  sign 
refers  to  quantities  on  the  transmitted  side  of  the  boundary.  Thus,  computations  with 
the  propagator  matrix  of  dynamic  ray  tracing  are  not  completely  analogous  to  those 
using  the  stress-displacement  propagator.  In  the  former  case,  the  dynamic  ray  tracing 
propagator  is  discontinuous  at  first  and  second  order  velocity  discontinuities,  while  in 
the  latter  case  the  stress-displacement  propagator  is  continuous  at  solid-solid 
boundaries.  Nevertheless,  the  multiplication  of  dynamic  ray  tracing  propagators  is  still 
a  useful  device  for  extending  solutions  for  the  P  and  Q  matrices  from  2-  and  3-D 
regions,  in  which  it  is  necessary  to  numerically  integrate  the  dynamic  equations 
through  thick  1-D  and  homogeneous  regions,  for  which  the  P  and  Q  elements  are  known 
analytically.  The  patching  of  the  solution  is  most  conveniently  performed  at  pseudo 
boundaries  at  which  velocity  and  its  first  spatial  derivatives  are  continuous.  Patching 
can  be  performed  at  more  general  boundaries  provided  that  one  first  uses  the  jump 
conditions  to  correct  P.Q  elements  to  values  appropriate  for  the  transmitted  side  of  the 
boundary. 

In  order  for  the  n  matrix  to  have  the  properties  of  a  propagator,  one  must  use  the 
plane  wave  initial  conditions  for  P,  Q  of  Cerven^  (19B5a)  rather  than  the  modified  plane 
wave  conditions  recommended  by  Madariaga  (1984)  for  a  source  in  a  region  having  a 
non-zero  velocity  gradient.  The  modified  initial  conditions  of  Madariaga  (1984)  were 
proposed  to  obtain  phase  fronts  at  the  source  that  were  equivalent  to  plane  waves.  This 
condition  was  necessary  to  stabilize  the  computation  of  Gaussian  beam  seismograms.  It 
also  makes  Gaussian  beam  synthesis  more  closely  resemble  the  process  of  plane  wave 


superposition  upon  which  Maslov  synthesis  (Chapman.  1982)  is  formulated.  Recently, 
however,  Cerveny'  (1985b)  has  proposed  optimal  beam  width  parameters  that  are 
equivalent  to  plane  wave  fronts  at  the  receiver.  These  parameters  are  optimal  in  the 
sense  of  minimizing  the  error  in  the  truncation  of  a  beam  sum.  This  choice  of  beam 
parameters  also  addresses  the  stability  problem  noted  by  Madariaga.  For  a  source  and 
receiver  both  at  the  surface.  Cerveny'’s  optimal  beam  width  parameters  are 
reciprocally  equivalent  to  Madariaga's  modified  plane  wave  initial  conditions.  This  can 
be  demonstrated  using  the  propagator  matrix  of  dynamic  ray  tracing  and  reciprocal 
relations  between  elements  of  the  P  and  Q  matrices. 

6.  ITERATIVE  PARAXIAL  RA  Y  TRACING 

As  an  additional  check  on  Gaussian  beam  synthetics,  it  is  planned  to  also  construct 
synthetics  using  the  KirchhofT  integral  technique  (Haddon  and  Buchen,  1981;  Scott  and 
Helmberger,  1983).  In  this  technique  an  integral  must  be  evaluated  on  a  reference 
surface  between  source  and  receiver.  For  heterogeneity  near  the  receiver  site,  this 
reference  surface  is  most  conveniently  selected  at  the  base  of  the  heterogeneity. 
Numerical  integration  of  the  dynamic  and  kinematic  ray  tracing  equations  can  be 
performed  in  a  model  of  the  structure  beneath  the  source  site.  Analytic  forms  can  be 
used  for  the  spreading  and  travel  times  needed  for  the  long  paths  in  a  radially 
symmetric  model  between  the  reference  surface  and  the  receiver.  It  is  planned  to 
perform  .ins  surface  integral  by  a  parameterization  of  the  surface  in  terms  of 
isochrons  (contours  of  equal  arrival  time  between  source  and  receiver,  as  described  in 
Haddon  and  Buchen,  1981).  An  intermediate  step  to  the  determination  of  isochrons  is 
to  densely  tabulate  travel  time  from  the  source  point  to  the  reference  sur  face.  To 
perform  this  it  is  planned  to  use  an  iterative  ray  tracing  technique.  In  this  scheme  the 
paraxial  approximation  is  used  to  calculate  take-off  angles  of  rays  connecting  the 
source  point  with  a  particular  point  on  the  reference  surface.  This  scheme  is  described 
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in  f^erveny  {19B5a),  and  amounts  to  estimating  the  new  angles  71,73  needed  to  hit.  a 


point  by  the  corrections 


where 


7ln«»u|  _  |7loJ<i|  ^ 

72  nnu  72  old  A72 


‘si  = isl 


with  (91, 9a)  being  vector  between  the  point  desired  and  the  point  actually  hit,  measured 
in  component  directions  of  the  local  ray  centered  co-ordinate  system.  With  this 
correction  it  is  possible  to  shoot  a  sparse  fan  of  rays  onto  the  reference  surface  and  to 
interpolate  travel  times,  spreading  functions,  and  take-off  angles  (needed  for  the 
source  radiation  pattern)  between  densely  spaced  grid  points.  An  example  of  the  rapid 
convergence  of  eq.  5  is  shown  in  Figure  4  for  a  3-D  model  and  initial  guesses  for  the 
take-off  angles. 


5.  CONCLUSIONS  AND  FUTURE  WORK 


A  test  structure,  determined  from  block  3-D  inversion  of  teleseismic  travel  times, 
gave  amplitude  variations  of  about  a  factor  of  2  along  a  100  km  long  profile.  The  block 
sizes  of  the  test  structure  varied  from  10-20  km.  in  linear  dimension  with  velocity- 
variations  on  the  order  of  2  to  8  per  cent. 

Future  work  will  be  directed  toward  comparing  the  amplitude  fluctuations 
predicted  by  3-D  models  having  varying  degrees  of  resolution  in  order  to  determine  the 
resolution  required  to  obtain  significant  amplitude  corrections.  A  parallel  series  of 
experiments  will  be  conducted  on  synthetic  velocity  models  constructed  to  have  varying 
distributions  of  characteristic  scale  lengths  and  velocity  fluctuations  within  the  domain 
of  validity  of  the  Gaussian  beam  method. 


Synthetics  will  be  constructed  at  teteseismic  distances  using  beam  superposition, 
ray  theory,  and  the  Kirchhoflf  integral  technique.  A  key  step  in  these  calculations  will 
be  the  patching  of  dynamic  solutions  between  1-D  and  3-D  portions  of  the  earth  model 
using  the  propagator  matrix  of  dynamic  ray  tracing.  This  will  enable  3-D  model 
variations  to  be  easily  included  at  both  the  source  and  receiver  ends  of  the  ray  paths. 
Thus  the  combined  effects  on  amplitudes  of  source  and  receiver  structures  can  be 


examined. 
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Figure  la.  Ray  plot  of  downward  propagating  rays  for  a  shallow 
focus  event  in  the  Zandt  (1981)  structure. 

Figure  lb.  Surface  project  of  the  rays  shov/n  in  Fig.  la. 

Figure  2a.  Rays  from  a  teleseismic  plane  wave  vertically  incident 
on  the  Zandt  (1981)  structure. 

Figure  2b.  Same  as  Fig.  2a;  denser  ray  spacing. 

Figure  3a.  Ray  intersections  with  the  surface  for  the  incident 
plane  wave  described  in  Fig.  2a  and  2b. 

Figure  3b.  Gaussian  beam  synthetic  seismograms  along  the  profile  AA' 
in  Fig.  3a. 

Figure  4.  An  example  of  the  convergence  of  iterative  paraxial 
ray  tracing.  Side-view  and  map-view  of  ray 
trajecto'-ies  are  shown.  In  the  side  projection  the  target 
points  are  at  -3.5,  -4.5,  and  -5.5  km.  along  the  right 
hand  margin.  All  targets  lie  along  the  (3.0,  0.0)  km  co-ordinates 
of  the  map-view.  Note  that  even  a  poor  initial  guess  for  the 
take  off  angles  in  the  map-view  quickly  converges  to  the  target 
point.  The  target  points  lie  within  a  3-D  low  velocity  zone, 
which  perturbs  a  piecewise  vertical  gradient  model. 


Figure 
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